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The speed and the versatility of today's computers open up new opportunities to simulate 
complex biological systems. Here we review a computational approach recently proposed 
by us to model large tumor cell populations and spheroids, and we put forward general 
considerations that apply to any fine-grained numerical model of tumors. We discuss 
ways to bypass computational limitations and discuss our incremental approach, where 
each step is validated by experimental observations on a quantitative basis. We present 
a few results on the growth of tumor cells in closed and open environments and of tumor 
spheroids. This study suggests new ways to explore the initial growth phase of solid 
tumors and to optimize anti-tumor treatments. 
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1. Introduction 

1.1. Tumor growth as a multi-scale biological process 

Biological systems span vast spatiotemporal scales, from the microscopic dynamics 
of atoms to the macroscopic dynamics of cell clusters. Information flows in both 
directions and determines the behavior of living matter and ultimately the normal 
physiology of organisms and the onset of pathologies such as tumors. Individual 
tumors are complex biological systems and, in spite of great therapeutic advances, 
many tumors still escape treatment and lead to death. Indeed, the malignancy and 
the response of tumors to therapy depends on their growth potential which in turn 
is determined by the ability of tumor cells to adapt to different environments, to 
compete with normal cells for both space and nutrients and to ignore molecular 
signals attempting to block cell cycling or to promote cell death.i^ Major efforts 
have been made by experimental scientists to highlight the molecular circuits, that 
are often altered as the consequence of genetic alterations, underlying tumor cell bi- 
ology and, on this basis, to develop novel therapeutic strategies. This has produced 
a huge body of knowledge which has deepened our understanding of the molecular 
details of tumor cell biology but often with little or no consequence for the clinical 
management of tumors. Tumor resistance to radiotherapy and chemotherapy is in 
fact the main reason for ineffective therapy, and this lack of responsiveness is only 
partly caused by genetic alterations of tumor cells. Non genetic factors, such as 
pH, ions other than and oxygen distributions, nutrient supply, proteins of the 
extracellular matrix, that constitute the tumor microenvironment, have also been 
found important. Indeed, microenvironmental changes can protect tumor cells from 
apoptosis, promote their survival an d th e progression of both epithelial tumors and 
non-solid hematologic malignancies. 12121 

Unfortunately, it is very difhcult to investigate the dynamics of microenvironmental 
changes in real tumors, as this would require non-invasive time-resolved technologies 
with appropriate spatial and temporal resolutions, and even if these technologies 
were available it would be extremely complicated to explore competing hypotheses 
of tumor microevolution, to detail and control the underlying mechanisms and to 
predict the impact of this knowledge on the clinical management of tumors. 
At least part of the complexity of the problem is a sheer consequence of tumor size: 
clinicians deal with the macroscopic properties of tumors, i.e. masses that may even- 
tually weigh a few kilograms, and thus with a number of cells that ranges between 
10® and 10^^, and that may grow for months or years, with a corresponding num- 
ber of cell cycles somewhere in the range between 100 and 10000. Moreover, at the 
microscopic level the malignant transformation of single cells is a multistep process 
that involves the modification of several molecular circuits which, in turn, modify 
the cell's behavior and the relationships between cells and the environment.^ In 
addition, the epigenetic and environmental factors, which include cell-cell interac- 
tions, also conspire with genetic information to make tumor growth a highly variable 
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process with very strong feedbacks. The highly nonlinear character of the cells' 
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internal molecular machinery, combined with the cell-cell and environmental inter- 
actions, with the large number of cells in a tumor, and with the extended tumor 
lifespan, make predictions based on the behavior of individual molecular circuits 
utterly haphazard. 

The availability of powerful computers has already helped bridge the gap between 
observations and predictions in many complex problems, and this suggests that in 
the future we shall be able to simulate the behavior of large cell populations ab 
initio, starting from individual molecular reactions in single cells and climbing the 
ladder of complexity up to the behavior of whole multicellular organisms. This re- 
quires an incremental approach, a ladder of increasingly complex models of cells, 
bearing in mind that each step must be validated by experimental observations. 

1.2. Multicell tumor spheroids: an in vitro cell model with 

intermediate complexity between real tumors and conventional 
tumor cell cultures 

Small tumor cell aggregates (volume < Imm^) may escape conventional treatment 
of solid tumors and, in time, may lead to recurrence of the primary pathology, of- 
ten with a different phcnotype (e.g., acquired resistance to chemotherapy, acquired 
ability to metastatize, etc.) pl^' cells can also grow up to masses of this size without 
the support of the vasculature, although recent work is challenging this traditional 



Unfortunately it is very difhcult to study these micromasses both in humans and 
in animal models, because their size is below the imaging limit of current tech- 
nologies and it is not possible to measure their biological parameters and obtain 
information to validate the results of numerical simulations and draw conclusions 
on their biological and clinical properties. Multicellular tumor spheroids represent 
a valid and effective experimental cell culture techniqu e cap able of preserving the 
three-dimensional topology of actual tumor cell clusters.'-'-'"' Indeed, it is the three- 
dimensional topology that determines many important biological features, like the 
expression of specific genes, a slowed-down diffusion of nutrients and waste, and 
also the expression of new phenotypes like the resistance to radiotherapy, and in 
fact multicell tumor spheroids display many inte resti ng biological properties that 
cannot be observed in monolayer cultures such asP^ 

• heterogeneous expression of membrane receptors (that regulate cell adhesion and 
metabolism and also may act as target for specific anti-tumor drugs); 

• production of an intercellular matrix (important for cell aggregation and for pen- 
etration of cells of the immune system); 

• heterogeneous distribution of nutrients and oxygen that lead to the formation of 
a necrotic core and to a gradient of cell proliferation; 

• appearence of resistance phenomena and/or heterogeneous response to antitumor 
therapies; 

• growth kinetics very similar to those observed in vivo. 
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Multicell tumor spheroids are thus intermediate between traditional ceU cuhures 
and tumors in vivo, and at the same time they are accessible to experimental mea- 
surements: they provide many data that can be used to test and validate multi-scale 
models of solid tumor growth in the prevascular phase. They are morphologically 
similar to small tumors below the detection threshold, and they share with them 
the lack of vascularization. For these reasons tumor spheroids are the perfect targets 
for mathematical and computational models of tumor growth. 



1.3. Tumor spheroids as targets for mathematical and 
computational models 

Several mathematical and computational models of solid tumors and of spheroids 
have been developed since the pioneering work of Burton,'^ and most of them have 
singled out some important feature of tumors, such as the diffusion gradients of pH, 
oxygen and nutrients, or more recently, the formation of new blood vessels and/or 
the migration of tumor cells out of the primary lesion. It is not the aim of this work 
to provide a comprehensive review of the ma ny mo dels developed so far, since this 
has been the subject of recent survey papers .E^^li^ We note, however, that: 

• most models (with few exceptions like, e.g., refs! '^"^l^^l^^ |), assume arbitrary units 
and outputs cannot be compared to actual experimental data: they describe only 
qualitative behaviors that somehow resemble those displayed by real avascular 
tumors. Due to the nonlinear character of the underlying biological processes, and 
hence of model equations, many different models can in principle produce outputs 
describing the same qualitative behaviors, and it is not possible to validate/falsify 
such models; 

• the appearance of widely different resistance phenomena to antiturnor therapies 
in similarly grown, isolated, tumor spheroids of the same cell typ^^^ indicates 
that random fluctuation phenomena play an all-important role in the growth 
kinetics of spheroids. It is well-known that the discrete events at the single-cell 
level do display some randomness (e.g., mitosis, cell death, partitioning of sub- 
cellular organelles and molecules between daughter cells at division, etc.) and one 
can pinpoint the source of large-scale variability on these fluctuations, as they are 
amplified and propagated by cell-cell and cell-environment interactions; 

• to the best of our knowledge, experimental evidenc e on the migration of labelled 
cells and microspheres into tumor spheroidJ^^QS _ that reveal the existence of 
convective cell motions in the spheroid structure - have no t ye t been fully ex- 
plained by mathematical models in spite of various attemptsp^ 

• in real tumor spheroids, cells communicate with other cells and the environ- 
ment. In particular, concentration gradients of several molecular species, facil- 
itated transport processes into and out of individual cells, and the mechanical 
forces that push and pull cells as they proliferate with repeated mitoses and then 
shrink after death, play an extremely important role in influencing the biochem- 
ical/biomechanical aspects of spheroid growth, and they depend in turn on the 
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structure of the extracellular space. These processes mix with complex nonlinear 
interactions between the biochemical and the biomechanical part; 

The absence of vascularization in multicell tumor spheroids and their near 
sphericity hide an internal complexity which is not easy to tame either with an- 
alytic models, or with numerical models based on rough simplifications of the bi- 
ological settings such as cellular automata or lattice-based models. New modeling 
approaches should therefore be developed, and these should attempt to reproduce 
the biological complexity shown by real spheroids such as the interaction of cells 
with their neighbors and with the environment. Nonlinearity and discreteness of 
processes at the single-cell level suggest that the search for "logical rules" with sim- 
ple grid-based models will not be sufficient, and that models that are capable of an 
equivalent level of complexity may be required. Indeed, our final goal is the simulation 
of the microscopic interactions among cells and their environment to understand 
how these could affect spheroid growth and response to treatment. 

1.4. Complexity issues in real-life modeling approaches of cell 
populations 

In principle, an all encompassing simulation could start with the atoms in the cell: 
using the methods of molecular dynamics it would then be possible to simulate a 
living being starting from atoms, molecules, and a description of the forces that 



We can clarify the complexity issue by considering the problem of exploding memory 
size. If we take a cell radius '~ 5 fim, then the cell volume is approximately 5 x 
-|^Q-i6 j-|^3^ g^j^jj |-]-^g corresponding cell mass is about 5 x 10~^*'g, and this means that 
a single cell corresponds to about lO^'^-lO^'* atoms. On the other hand, if we take a 
spatial resolution 0.01 nm (approximately one tenth of the diameter of a hydrogen 
atom), and aim to simulate a system size of 1 mm, then for each coordinate we 
need at least a 24-bit dynamic range (3 bytes per coordinate), and thus at least 19 
bytes/atom (3 coordinates plus 3 velocities plus one atom label), and about lO^''- 
10^^ bytes/cell. Finally, the full simulation of a 10^ cell spheroid would require at 
least 10^°-10^^ bytes/spheroid, and we see that present-day computers are pitifully 
inadequate for such a brute force approach. 

Likewise, the time complexity of simulation algorithms would also be unmanageable: 
here we must assume that we are somehow able to tame the 0{N^) complexity 
of binary interactions between the simulation elements and also the complexities 
of several subalgorithms like matrix inversion and the like, and that the overall 
algorithmic time-complexity of the simulation program adds up to a mere 0{N). 
Since the fastest dynamics in acqueous solutions is determined by the motion of 
protons in the hydrogen bonds in water,!^ and is of the order of 1 ps, and since 
there are approximately 10^"^ hydrogen bonds in each cell, then one must take at 
least 10^^ time steps just to simulate 1 s of a single proton motion in a cell (and 
with a rather poor time resolution). 
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These approximate calculations amount to an operational definition of biological 
complexity, and they show that at the moment we cannot even dream of carrying 
out a true ab initio simulation of tumor spheroids. Therefore - in our approach ~ we 
have chosen the highest level scale at which the complexity of individual molecular 
processes is still somehow evident, the cell scale, although the question remains 
whether it is possible to stay at this level and still provide a comprehensive and 
realistic description of cell behaviour. 

2. Models 

2.1. A minimal model of the tumor cell 

When viewed at the mesoscopic scale the intracellular molecular machinery displays 
an astonishing complexity, with a huge number of intertwined chemical reactions 
that mark the different phases of the cell's life. On the other hand, a biophysical 
simulator of tumor spheroids must start from a realistic description of the tumor 
cell, and this ultimately means that at least cell metabolism and its interconnection 
with the cell cycle must be modeled at a sufficient level of detail in order to describe 
how the behavior of a single cell is affected by the other cells in the cluster, by 
the chemical composition of the environment and by physical parameters such as 
temperature, density, radiation, etc. 

Our approach is based on the fact that biochemical networks in the cell possess 
a hierarchical structurc^^HI jg known that, if a network has such topology, then 
the system dynamics are dominated by the network's hubs.^^Ol Thus, by modeling 
the hubs of the cell's biochemical network one should, at least in princip le, be 
able to capture most of the information of the cell's biochemical dynamics .EHH 
This also means that we must necessarily parameterize and average many details 
of cell metabolism. In this way, we achieve a significant reduction in computational 
complexity and a considerable reduction of the space-time scale problems that affect 
simulations aimed at calculating the properties of macroscopic objects starting from 
microscopic models. 

Fi gure 1 shows a sketch of the metabolic pathways that have been modeled 
so farPIESl •\Ye have taken into consideration both cell cycle-independent and cell 
cycle-dependent pathways. Among the former are the uptake and utilization of nu- 
trients such as glucose and glutamine and the elimination of waste molecules such 
as lactate. This part describes quantitatively the production of ATP by both oxyda- 
tive phosphorylation and glycolysis and it bridges the two paths: in f act, ATP and 
glutamine constitute the building blocks of proteins and DNA ^u^^ g^j^f^ their 
availability drives the kinetics of protein and DNA synthesis. These in turn deter- 
mine the expression of proteins that regulate the cell cycle such as pRb and cyclins 
and the duration of specific phases of the cell cycle such as the S phase .I^Sl Proteins 
that regulate the cell cycle also determine the duration of the other phases, and we 
have implemented a thresholding mechanism based on mu lt isit e phosphorylation to 
determine when a cell steps beyond cell cycle checkpoints. ^^^^^ 
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Fig. 1. Metabolic network implemented in the simulation program. See text for details. A complete 
description of this network can be found in reference 22 

ATP availability also determines the proliferation of mitochondria, a biological pro- 
cess that has been coupled phenomenologically to the increase of the cell volume, 
while ATP deprivation leads to cell death - although the model also takes into 
account several other mechanisms of cell death. Finally, we take into consider- 
ation stochastic aspects such as the random partitioning of mitochondria and the 
random distribution of molecules at cell divisi on, th at contribute to desynchronize 
the duration of the cell cycle in daughter cells .1^2^ 

The rationale for sorting the above processes out of the complex met abolic machin- 
ery of the cell has been discussed in detail in our previous papers .1^1^^ Because 
of the interconnections among the pathways, we believe that this model cannot be 
further reduced, above all if we are to produce a realistic simulation of how cells 
exchange biochemical information with the environment and the other cells, and 
regulate their life cycle. 

2.2. Modeling tumor spheroids 

The metabolic model of tumor cells shown in figure 1 is the basis on which we build 
the simulation of tumor spheroids and, for validation purposes with experimental 
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data, it has been used to sinrulate the growth of large ceh populations in both 
closed an d open environmentsl^i^^ The simulations of dispersed cells have been 
successfuljI^lEll however the modeling oftumor spheroids requires several other 
biological, chemical and physical aspects PSEOl 



(i) adequate numerical modeling of diffusion processes. The diffusion of chemicals 
in a cell cluster proceeds either by normal diffusion or by facilitated diffusion 
across cell membranes. Facilitated diffusion is mostly a biochemical process 
that has a weaker dependence on concentration gradients and is brought about 
essentially by transporters expressed at the cell surface. For facilitated diffusion 
to proceed in cell clusters, it is necessary to include the extracellular spaces 
as basic elements. Neighboring cells modify the surrounding environment and 
compete for the same resources, and thus it is important to define the proximity 
relationships among cells in the cluster. This is accomplished by a specialized 
part of the program (see below); 

(ii) environment. The environment is also taken into account in the simulation: it is 
defined as the external volume of nutrient fluid that communicates by diffusion 
with the extracellular spaces that surround cells, and it is modeled as a single 
compartment. The environment is modified both by the fast diffusion processes 
that transport nutrients and metabolites into and out of the cell cluster, and, 
in case, by the slow flushing of nutrient fluid in a bioreactor enclosure; 

(iii) biomechanical interaction. In a cluster, cells interact mechanically as well as bio- 
chemically; this part of the simulation program is essentially a simple integrator 
like those found in dissipative dynamics simulations. Ce lls a re approximated by 
soft spheres that move in a highly viscous environment 

(iv) geometry. Both the mechanical and the diffusion part of the program require 
a knowledge of the proximity relationships among neighboring cells. We wish 
to point out that, unlike other models, in our simulation program cells do 
not grow in an environment defined by a fixed grid, where it is difficult to 
model appropriately processes such as cell division, and that cells are free to 
move in the three-dimensional space as as they are pushed and pulled by forces 
resulting from biomechanical interactions between cells. The nearest neighbors 
are defined by the links in a Delaunay triangulatiort^ and they are computed by 
the triangulation methods in the computational geometry package CGAL (see 



http://www.cgal.org). In this way all the computational complexity of binary 
interactions is reduced from a potential O(iV^) to a much more manageable 
OiN). 



A full description of our model that includes all technical details can be found 
in reference 32. 
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2.3. Stability issues in the numerical integration of model 
differential equations 

In our model each cell is described by a metabolic network and by other mechanisms 
that include both discrete deterministic and stochastic events. The description is 
thus mixed, with smooth evolutions interspersed with discrete steps. The exchange 
of molecules with the surrounding environment means that transport into and out 
of cells is closely linked with diffusion processes that involve the whole cluster of 
cells, and finally lead to a very large set of differential equations. Since our goal is 
to simulate tumor spheroids up to a diameter of 1 mm, which corresponds to about 
1 million cells, the software must eventually solve a very large number of coupled 
nonlinear equations, as many as 10''-10^ equations (because there are - at present 
- 25 variables per cell) . These equations are similar to other equations encountered 
in systems biology, but the number is uncommonly large. In addition, the equations 
are quite stiff, since they describe processes that range from fast diffusion in small 
extracellular spaces (a few tens of /is) up to the slow development of the spheroid as 
a whole 10^ s) and thus the characteristic times span about 12 orders of magni- 
tude. We have solved the complex stability problems that arise in such a situatiorP^l 
and thus this model now stands up as a true multiscale model, both in space and, 
even more so, in time. We argue here that a modeling such as that described in 
reference 33 is mandatory in any such simulation of cell clusters. 
We also remark that some of the model parameters slowly change as cells grow and 
this is once again at variance with most differential systems used in systems biol- 
ogy where parameters are fixed. Finally, the continuous system evolution described 
by the differential system is interrupted at random times by discrete events; these 
events may be internal transition in individual cells (in this case the system pa- 
rameters change abruptly from one integration step to the next) or cellular division 
events (in this case the number of equations changes). For all these reasons, we 
believe that our simulation program is unique in the sense that it tackles simultane- 
ously and for the first time a vast array of technical issues that are not addressed by 
any other model. In particular, the same implementation scheme of model equations 
can be used to simulate both fast and slow processes simply by tuning the integra- 
tion time step, and this is an added value to the simulation program that might be 
exploited to investigate at will biological events with different characteristic times. 

3. Results 

3.1. Simulations of the growth of dispersed cells in a closed 
environment 

One important aspect of the approach described in the previous section is that it 
is fully quantitative and its outputs can be directly compared with available ex- 
perimental data. The model of cell metabolism, growth and proliferation described 
above can be used to simulate a population of dispersed cells growing in a closed 
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environment. This is equivalent to considering cultures of blood cells in vitro, such 
as leukemia cells. Cells growing in a closed environment establish a sort of nega- 
tive feedback with the environment itself. While cells grow, they consume nutrients 
and release waste molecules that acidify the medium. As the environment gradu- 
ally becomes more and more acidic, the uptake of nutrients is also reduced and 
can eventually switch off completely, thereby leading to a depletion of the energy 
reserves and ultimately to cell death. This mechanism involves the whole model 
of cell metabolism and control of the cell cycle and can be tested experimentally 
because it defines the carrying capacity of the environment where cells are grown. 
We have already demonstrated that our model nicely reproduces common growth 
curves observed in vitrH^ and shows predictive capabilities, in that it can also de- 
scribe the growth of leukemia cells under non conventional conditions, where, e.g., 
the biochemical composition of culture media is periodically refreshed. The model 
predictivity in nonconventional conditions is obviously an important test, because 
it demonstrates that the model is not just a qualitative description, but can be used 
to set up a true in silico laboratory. Indeed, the model outputs are in good quan- 
titative agreement with actual data on metabolic and cellular parameters (Table 



3.2. Simulations of the growth of dispersed cells in an open 
environment 

The negative feedback between cells and their environment discussed above can be 
partially removed by opening up the environment, setting up a flow that removes 
consumed culture medium and replaces it with fresh medium. Experimentally, this 
condition is realized in bioreactors such as those used to culture cells at high density 

OA 

for biotechnology purposes, e.g., antibody production.- Under these conditions, vi- 
able cells are expected to reach a steady state given by the dynamic equilibrium 
between proliferation and death. 

Figure 2 shows the result of a simulation campaign carried out with the numer- 
ical model in reference 22. Parameter values were left unchanged with respect to 
previous simulations, and the only difference is that here a continuous flux of fresh 
medium slowly replaces the culture medium with rates comparable to those in actual 
experiments with bioreactors. As we see in figure 2, there is a good correspondence 
between simulations and actual data for different flow rates. 

This kind of virtual experiment was not planned during model development and 
model parameters were not tuned to take into account the growth of simulated 
cells in an open environment. This further demonstrates the predictive power of the 
model and its potential use as an in silico cell biology laboratory. 
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Table 1. Estimated morphologic, kinetic and metabolic parameters for a population 
of dispersed tumor cells and comparison with actual experimental data. 



Parameter Simulated Experimental Reference 



Morphologic 


J^verage 




Max 






38, 


Radius {fim) 


5.0 


4.8 


5.3 


5.5 


- 7.1 


Volume(/im^) 


530 


471 


623 


700 


- 1500 


,38.391 


Mitochondria/cell 


220.4 


190.6 


266.9 


83 


- 677° 


So] 


Kinetic 














Growth rate''(h-l) 


0.035 






0.03 


- 0.035"= 


122 


Doubling time^(h) 


19.8 






19.7 


- 22.8"= 


121 


Gl (%) 


52.5 


48.4 


59.3 


54.4 


± 2.2"= 


121] 


S {%) 


34.5 


30.5 


40.5 


27.5 


± 5.8"= 


[21] 


G2/M (%) 


12.9 


7.3 


17.7 


16.4 


± 1.7^ 


121] 


Metabolic 












121] 


ATP/cell'^ 


5.5 


5.4 


5.6 


4.3 


- 5.8 


Glucose uptake^ 


1.9 ± 0.3 






2.5 


± 0.2 


HI] 


Lactate production^ 


3.8 ± 0.3 






3.9 


± 0.8 


HI] 


ATP production'^'-'' 


19.8 ± 8.3 






37.8 


m] 


ATP production'^'f 


10.6 ± 1.3 






11.4 


± 2.3 


m] 


Oxygen consumption^ 


0.25 ± 0.1 






0.48 


± 0.1 





° Range of the number of mitochondria observed in different cell types 
''The growth rate for both simulated and experimental cell populations was calculated 
by exponential fitting of growth curves. The doubling time was then calculated as 
log 2 / (growth rate) 

^Data measured for MOLT3 (human T lymphoblastoid cell line) and Raji (human B 
lymphoblastoid cell Hne) cells in our own experiments 

Values are expressed as 10 kg 
^Values are expressed as 10~^® kg s~^ 
^ATP production through oxidative phosphorylation 
^ATP production through glycolysis 

3.3. Climbing the third dimension: a simulation test run of the 
growth of tumor spheroids 

Our model can be used to simulate the behavior of large tumor cell populations 
growing both in closed or in open environments. The very same model can also 
be used to simulate the growth of avascularized tumor spheroids, where we have 
to take into account additional biological, chemical, physical and mathematical 
aspects. Cells in a spheroid adhere to one another. Thus, one must include and 
model biomechanical forces that act upon cell s. Cel ls have been modeled as soft 
spheres interacting through visco-elastic forces .122122] ^pj^ege forces allow the whole 
cell aggregate to preserve a three-dimensional structure in spite of the push that 
cells exert at mitosis on neighboring cells as they compete for space. This results in 
a dynamic balance between repulsive and adhesion forces that shapes the spheroid 
structure in the course of its development. 
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20 5 10 15 



Time (days) 

Fig. 2. Growth kinetics of dispersed cells in a bioreactor. Points represent experimental data taken 
at different flux rates of fresh medium and are: squares, 2.33 volumes of fresh medium perfused 
/ effective suspension volume / day (vvd); circles, 1 vvd; triangles, 0.48 vvd. Lines represent 
simulation outputs obtained with our program for the same fluxes. See also the text for further 
details. Experimental data are taken from reference 34. 

Each cell is surrounded by a small free volume that corresponds to the extracellular 
space. This is fundamental to allow nutrients and waste molecules to diffuse freely 
through the cell cluster. The inclusion of diffusion, however, introduces processes 
that occur with very short characteristic times, in the order of a few tens of /iS, and 
this increases the stiffness of the underlying system of differential equations. 

Figure 3 shows a result obtained in a run with a simulated spheroid of more than 
200000 cells and 16 days old (simulated time). The oxygen concentration decreases 
towards the center of the spheroid and one can easily detect a spatial gradient that 
is very similar to that observed in real experiments. As the spheroid grows the in- 
ternal environment become hypoxic, and this contributes to cell death and to the 
formation of a necrotic core. 

Once again, we stress that our approach is quantitative and this allows us to test 
model outputs with actual experimental data. The profiles of the main metabolic 
actors have been measured in spheroids and data on at least some cell and mor- 
phologic parameters are also available for comparison. In Table 2 we show that our 
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Fig. 3. Oxygen concentration profile in a simulated spheroid after 16 days of simulated time. 
Panel a; slice cut through the center of the spheroid showing the oxygen concentration. The gray 
levels map the concentration value: white, higher concentration; black, lower concentration (the 
background is light gray to provide suitable contrast). Panel b: oxygen concentration given as 
atmospheric partial pressure vs. the distance r from the centroid of the simulated tumor spheroid 
(the rightmost value corresponds to the standard atmospheric partial pressure). 

niodel of multicell tumor spheroids provides outputs that compare favorably with 
available microscopic data, and in Figure 4 we show that the growth kinetics of 
spheroids can also be properly simulated. 

Thus, our model connects the mesoscopic and the macroscopic scales. It is also 
important to note that the data in Figure 3, Figure 4 and Table 1 do not consti- 
tute a final validation of our model and more virtual experiments must be carried 
out, possibly under conditions that were not considered during model development. 
These include the growth of spheroids in media with a different concentration of 
glucose and oxygen. Experimental data are availabld^, and virtual experiments 
can be carried out because of the versatile model of cell metabolism that we have 
implemented in our model of multicell tumor spheroids. 

Since concentration profiles are computed during the time evolution of the whole 

cluster, the dynamic variation of the distributions of molecular species can also 

be appreciated and studied for the first time. This computational approach could 

eventually help to improve our knowledge on the initial growth phase of avascular 

solid tumors, as it discloses details that are otherwise very difficult to detect and 
(321 

measure. '^''^ 

4. Discussion and Conclusions 
4.1. Parameter count 

We model avascular tumor spheroids starting from individual molecular reactions in 
single cells and we climb the ladder of complexity up to large multicellular clusters. 
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Table 2. Estimated metabolic, histologic and kinetic parameters for a virtual 
multicell tumor spheroid and comparison with actual experimental data. 



Parameter 


Simulation 


Experiments 


References 


Metabolic 








Glucose uptake°(kg s"^ rD~^) 


1.44 ■ 10-3 


5.4 - 12.6 ■ 10-3 


ES] 
US] 


Lactate release''(kg s" m" 


1.35 ■ 10-3 


5.4 - 9 • 10-3 


p02''(mmHg) 


7 


0-20 


|43] 


pH^ 


6.7 


6.6 - 6.99 


|44|45| 


ApH'^ 


0.77 


0.49 ± 0.08 


HH 


Histologic 






I45I46I4VI 


Viable cell rim thickness^ (/im) 


155 


142 - 310 


Hypoxic rim thickness-^ (/^m) 


98 


44 ± 52 


Ha 


Kinetic 








Cell cycle distribution 








Gl (%) 


57.3 


58 ± 4 


US] 


S (%) 


21.6 


19 ± 1 


BS] 
is] 


G2/M (%) 


21.1 


23 ± 1 



Note: Metabolic and histologic parameters in spheroids of approximately 500 
^m diameter. 

° Rate of glucose uptake or lactate release per viable spheroid volume 

Central p02 tension (experiments) or estimated in the centroid (simulations) 
'^pH has been determined in the central region of the spheroids. This corre- 
sponds to a sphere radius ~ 100 ^m about the centroid of the spheroid 

Difference between environmental pH and pH 200 ^m below the spheroid sur- 
face 

"^In our simulations the viable cell rim thickness corresponds to the distance 
between the spheroid surface and the inner shell where only 5 % of the cells are 
still alive 

These values corresponds to the radius of the necrotic core 



As shown above, this imphes that we must deal with very different space and time 
scales, with the latter spanning 12 orders of magnitude. This effort unavoidably 
translates into a complex mathematical description that involves a vast number 
of equations and parameters, and with many free parameters one can in principle 
simulate any number of different patterns. Thus, a central question is whether our 
model can indeed be considered a realistic description of metabolism and mechani- 
cal evolution of cell clusters. On the whole the present simulation is defined by 100 
parameter values, and a realistic simulation requires a reliable set of parameters: 
to this end we extensively searched the scientific literature to find experimental 
measurements for as many parameters as possible. Whenever experimental mea- 
surements are missing, we assume values estimated from independent biophysical 
modeling of experimental data. Once fixed, the parameter values are not changed 
any more and simulation outputs are compared with new sets of experimental data 
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Time (hours) 



Fig. 4. Growth kinetics of a simulated muticell tumor spheroid and comparison with experimental 
data. The simulation (black line) was started by placing one virtual cell in an environment of 1 
ml filled with standard medium. The initial concentrations of glucose, glutamine, oxygen and the 
pH were set at the same values as for real experiments with tumor cells grown in a C02 incubator 
and in RPMI medium supplemented with 10 % fetal bovine serum. The data in the figure (black 
symbols) are the mean ± standard deviation radii of nine spheroids obtained with 9L cells (rat 
glioblastoma cell line). Spheroids were grown isolated in individual culture wells under standard 
conditions, and their size was measured daily under the microscope (this experimental set was 
described in reference 42). 



to test the predictive power of the model itself (as we have shown here for cells grow- 
ing in a bioreactor) . All comparisons between simulation outputs and experimental 
data are finally carried out on a strict quantitative basis. Thus we feel confident that 
our numerical simulator is indeed a reliable model of the growth of tumor cells and 
tumor cell clusters. Interestingly, we found that most parameters are correlated and 
can by no means assume arbitrary values. This is mainly due to the strong feed- 
back between cells and the environment where simulated cells grow and that we 
have modeled. For example, one might tune metabolic parameters to allow cells to 
utilize nutrients more efficiently, produce energy under the form of ATP and grow 
faster. But this consequently results in a higher production and secretion of waste 
molecules that increase rapidly the acidity of the medium thus leading to cell death. 
Many parameters of the metabolic network of the cell are also strongly correlated, 
and this is a typical consequence of the interconnections between reactions that 
utilize different substrates for the same purpose, such as in the case of ATP and 
glutamine for both protein and DNA synthesis. As a consequence, the actual size 
of the parameter space is greatly reduced, and parameter tuning is far less complex 
than it appears to be at first sight. 
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4.2. A real-life simulation 

Here we have reviewed the development of a numerical tool to simulate realistically 
the growth of avascular tumors and thus to explore the initial growth phase of solid 
tumors. This kind of numerical simulation has several important implications: 

(i) it is possible to perform virtual experiments in silico that complement in vitro 
measurements, where many parameters are not directly accessible, and also in 
vivo observations, where accessibility problems are even greater also because of 
ethical issues. Our simulation program is indeed a virtual laboratory where one 
can make experiments at will in due time, and we hope that in the near future 
it will drive experimentalists towards the search of yet unexplored biological 
properties of tumors; 

(ii) the simulation focuses the modeling effort on the important details of cellular 
biophysics and spawns new ideas, both theoretical and experimental. For exam- 
ple, one important aspect that we have investigated to test the validity of our 
simulation program is whether, and eventually under which conditions, it could 
simulate a cell population with desynchronized cell cycles as it is observed for 
real cells.^^H^ This prompted us to consider the relevant sources of internal 
randomness in cells, and in the attempt to investigate the biological causes of 
cell cycle desynchronizatio n we d eveloped theoretical tools and carried out new 
experimental observationsj^SEZl 

(iii) the numerical model includes many complex non-linear interactions between 
different parts of the cell, and thus it has interesting predictive properties as 
unexpected biological behaviors can emerge; 

(iv) the model integrates several parts of our knowledge of cell biology, a knowledge 
that is fragmented in a huge number of small pieces throughout the scientific 
literature. An important aspect of our effort is that we are trying to connect 
together at least part of these pieces, and check their overall consistence. The 
model produces results that compare favorably with experimental data, and this 
indicates that it is possible to understand the cells' functions at the systemic 
level in quantitative terms; 

(v) because of its incremental structure, our simulation program may serve as a 
platform to test the validity of other models of specific biochemical circuits. 

For these reasons, we believe that our effort is by no means just a modeling 
exercise, but a serious and novel attempt to model cell biology, and a bridge be- 
tween the molecular level - albeit often phenomenological - and the macroscopic 
observations, that demonstrates the possibility of carrying out realistic simulations 
and the evaluation of macroscopic parameters from microscopic descriptions. 
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